template<class Resample>
class Maruyama{
	public:
	Maruyama(Resample *R, mat X, mat Y)
	{
		_R=R;
		_p=X.n_cols;
		_n=X.n_rows;
		mat m(_p,1);
		mat m2(_n,1);
		m.fill(0);
		m2.fill(0);
		mat S(_p,_p);
		mat S2(_n,_n);
		S.eye();
		_X=X;
		_Y=Y;
		S2.eye();
		mat I;
		I.eye(_n,_n);
		_Xy=(2*diagmat(Y.col(0))-I)*X;
		Psy=S2-_Xy*inv(_Xy.t()*_Xy)*_Xy.t();
		_G=new Distribution::Gaussian(_p,m,S);
		_G2=new Distribution::Gaussian(_n,m2,S2);
	}
	~Maruyama()
	{
		delete _G;
	}
	mat operator()(int M)
	{
		mat Z=(*_G2).scrambled(M);	
//		mat Z3=(*_G2).scrambled(M);	
	//	cout << Z; 
		mat h(M,_n);
		mat b(M,_p);
		mat t(M,1);
		mat v(M,1);
	
		for(int i=0;i<M;i++)
		{
		//	t(i,0)=dot(Z3.row(i),Z3.row(i));
			h.row(i)=abs(Z.row(i))/sqrt((dot(Z.row(i),Z.row(i))));
			mat temp=h.row(i)*Psy.t();
			v(i,0)=(double)1/sum(sum(abs(temp)));
		}
		double *w=new double[M];
		double max=_n*log(v(0,0));
		for(int j=0;j<M;j++){
			 
			w[j]=_n*log(v(j,0));
			if(w[j]>max)
			{
				max=w[j];
			}	
		       //	cout << " " << w[j];
		}
		double sum=0;
		for(int j=0;j<M;j++){
			w[j]=exp(w[j]-max);
			sum+=w[j];
		}
		double sum2=0;
		for(int j=0;j<M;j++)
		{
			sum2+=pow(w[j]/sum,2);
		}
		cout << "ess " << (double)1/sum2;
		mat udt=growingvect(M);
		(*_R)(&M,w,udt);
		//cout << udt;
		Arangemat(h,udt);
		//Arangemat(v,udt);
		for(int i=0;i<M;i++)
		{
			mat temp2=h.row(i)*Psy.t();
			boost::random::gamma_distribution<> Gamma((double)_n/2,2/dot(temp2,temp2));
			rgamma= new RandomG::Random<boost::random::gamma_distribution<> >(Gamma);
			t(i,0)=(*rgamma)();

		}
		delete[] w;
		mat Z2=(*_G).scrambled(M);	
		mat XX=inv(_X.t()*_X)*_Xy.t();
		mat XC=chol(inv((_X.t()*_X)));
		for(int i=0;i<M;i++)
		{
			mat t1=h.row(i);
			mat t2=Z2.row(i);
			mat t3=sqrt(t(i,0))*XX*t1.t()+XC.t()*t2.t();
			b.row(i)=t3.t();
		}
		return b;


	}
	private:
	Resample *_R;
	int _p;
	int _n;
	mat _Xy;
	mat Psy;
	mat _X;
	RandomG::Random<boost::random::gamma_distribution<> > *rgamma;
	mat _Y;
	Distribution::Gaussian *_G;
	Distribution::Gaussian *_G2;


};

